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GRAVITATIONAL COLLAPSE OF A TURBULENT VORTEX 

WITH APPLICATION TO STAR FORMATION 

by Robert G. Delssler 
Lewis Research Center 
Cleveland, Ohio 44135 

ABSTRACT 

The gravitational collapse of a rotating cloud or vortex is analyzed 
by expanding the dependent variables in the equations of motion in two- 
dimensional Taylor series in the space variables. It is shown that the 
gravitation and rotation terms in the equations are of first order in the 
space variables, the pressure gradient terms are of second order, aud 
the turbulent viscosity term is of third order. The presence of a tur- 
bulent viscosity insures that the initial rotation is solid -body- like near the 
origin. The effect of pressure on the collapse process is found to depend 
on the shape of the initial density disturbance at the origin. Dimension- 
less collapse times, as well as the evolution of density and velocity, 
are calculated by solving numerically the system of nonlinear ordinary 
differential equations resulting from the series expansions. The axial 
inflow plays an important role and allows collapse to occur even when 
the rotation is large. An approximate solution of the governing partial 
differential equations is also given, in order to study the spacial distri- 
butions of the density and velocity. 
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I. INTRODUCTION 

Gravitational instabilities appear to play a dominant role in the for- 
mation of stars and other astronomical objects, and calculations of the 
gravitational collapse of clouds (protostars) have been carried out by 
a number of authors (e. g. Larson 1969, 1972, 1973; Disney et al. 1969; 
Pension 1971; Tscharnuter 1975). Most of these have been numerical 
solutions of the governing partial differential equations, where, except 
in recent work such as that of Larson (1972) and of Tscharnuter (1975), 
the effects of rotation have been neglected. But as will be seen, the 
effects of rotation are of the same order as those of gravity, and are of 
lower order than those of pressure gradients. In the past work boundary 
conditions were assumed at a hypothetical outer boundary, where con- 
ditions were generally not well known. 

The present treatment differs from previous work in that the depen- 
dent variables are expanded about the origin in truncated power series 
in the space variables. This converts the governing partial differential 
equations to ordinary differential equations in time. Aside from the fact 
that the resulting ordinary differential equations are much easier to solve 
than are the original partial differential equations, this procedure has 
the advantage that the various physical processes are conveniently sepa- 
rated into first, second, and third order effects. Moreover the introduc- 
tion of boundary conditions at an outer boundary is replaced by the natural 
assumption that the dependent variables and their lower order derivatives 
are finite at the origin. That is, the boundary conditions are applied at 
the origin, and we do not have to specify the extent of the gaseous cloud. 
However, we still have to specify the size of the initial disturbance. 
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The effects of pressure are included and are found to depend on the 
amplitude and shape of the initial disturbance. The effects of rotation 
and of turbulent viscosity are also included in the analysis. Although 
the quantitative effects of turbulent viscosity are only of third order, it 
appears that the presence of the turbulent viscosity is important for in> 
suring that the velocity and its derivatives are bounded at the center of 
contraction and rotation. 

Most of the calculations were carried out by using ordinary differen 
tial equations in the independent variable time, as described above. 
However, in some of the calculations the radius was also retained as a 
variable in order to investigate the spacial variations of density and 
velocity. 

n. BASIC EQUATIONS 

The equations of motion and continuity for an axially symmetric 
compressible flow can be written in cylindrical coordinates r, 6 and 
z as 


0t r dr dz dr p dr 

(1) 

~ - w-^ - D 
9t 3r r dz 

(2) 

~ ^ 3w d<p 1 dp 

9t 9z 9r 9z p 9z 

(3) 

and 


= . 1 i. (rpu) - i- (pw) 
9t r 9r 9z 

(41 


where the gravitational potential tp is given by the Poisson equation 
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=1 J.4iS?V 4ffGp (5) 

az2 

and where u, v, and w are the velocity components in the r, 0 and 
z directions respectively, t is the time, p is the density, p is the 
pressure, G is the gravitational constant, and D is a turbulent or vis- 
cous drag term in the 0 direction. Drag terms are not shown in equations 
(1) and (3) because those terms have been assumed small compared with 
(or included in) the pressure gradient terms. To relate the pressure 
to the density, we use the polytropic relation 

P = PjPi'^P^ (6) 

or 

= rPiPi'^p*^"^-^ (6a) 

3r ^ ^ 3r 

where y (assumed constant) is the polytropic exponent for the collapse 
process. The subscripts one designate ambient values. 

The drag term D arises mainly because of the effects of turbulence, 
the effect of molecular viscosity usually being comparatively small for 
the high Reynolds numbers in astronomical systems. For our present 
purposes it should be sufficiently accurate to represent the effects of 
turbulence by a uniform turbulent viscosity c as in Deissler and 
Perlmutter (1960). Thus we write 



In Deissler and Perlmutter it is supposed that € is determined by the 
shear, and an estimate of its value is given by using a modification of 
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von Karman^s similarity theory. It is shown there that for v propor 
tional to r”^ (large radial flow) 



^ 1^1 


( 8 ) 


where k is the Karman constant, and r^ and are respectively the 
radius and tangential velocity at the outer edge of an initial disturbance 
to be specified. From the experiments cited in Deissler and Perlmutter, 

/c^/2»l/20 (9) 

In the present case, V lies between values given by an r'^ and an 
r variation, so that c as determined by shear will tend to be lower than 
the value given by equation (8), We will retain equation (8) in the present 
study as an upper limit for £ as determined by shear. Other effects such 
as normal strain (Deissler 1968, 1972) and gravitational instabilities will 
tend to offset the decrease in c associated with decreased shear, so that 
equation (8) may give a reasonable estimate. As will be shown, the 
effects of turbulent viscosity are quite small so that the exact value used 
for € is not critical. 

The set of equations (1) to (9) is determinate, and its solution will 
be considered in the next section. 

in. SOLUTION BY TAYLOR SERIES 
We can expand the dependent variables, u, v, w, p and (p (repre- 
sented by x) in two-dimensional Taylor series about r = z = 0, truncated 
after terms of third order in r and z, as 

X = Xo + Xj.r 4 X55Z + i 4 x^^rz 4 1 


6 


1 

6 




.3 . J. 
2 




(10) 


where 


Xq = (x)i.,2> ’‘r = f— > \z =f— 1 • 

-n [dr/r=z \BrSz/r=z 


etc. 


=0 =0 
We take the r and z axes as axes of symmetry, so that the dependent 
variables are symmetric or antisymmetric about those axes. Then 
(see Fig. 1) 


w(r, -z) = -w(r, z) 
w(-r, z) = w(r,z) 
w(-r, -z) = -w(r, z) 
<p(r, -z) = (p(r, z) 
<;p(-r,z) = <p(r,z) 
cp = (-r, -z) = (p{r,z) 


u(r, -z) = u(r, z) 
u(-r, z) = -u(r, z) 
u(-r, -z) = -u(r, z) 
p(r, -z) = p(r, z) 
p(-r, z) p(r, z) 
p(-r, -z) = p(r, z) 

The vortex rotates about z, so that 
v(-r,z) = -v(r,z) 
v(r, -z) = v(r, z) 
v(-r, -z) = -v(r, z) 

If we impose these symmetry conditions on equation (10), we get 

u = u„r + -i u^^^ r^ + -i u^„^ rz^ 
r 0 rrr 2 

V = v^r + i v^^^ r^ + -i v^^„ rz^ 

r 0 rrr 2 

« « . ^ 1 « ^ 2 , ^ 1 ^ ,3 

w = w_z + — w^^_ r z + — w___ z 

z 2 g 


( 11 ) 

( 12 ) 

(13) 




(14) 
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(15) 

Consider first the gravitational potential <p. Substituting equa- 
tions (14) and (15) in (5) gives, for coefficients of r^ 

2<?>rr + ’’zz = 

The relation between and can be considered as a boundary 
condition for (p. For that boundary condition we set - (p^^; that is 
we assume spherical symmetry for cp near the origin. Larson (1972) 
found that his results were not sensitive to the particular boundary condi- 
tions used for (p. Thus equation (16) becomes 

^rr “ ^zz ~ 

From equations (15) and (17), 


and 


-^ = -Gp„r 
3r 3 " 


= 4i Gp.z 

0z 3 


(18) 


(19) 


Substituting equations (11) to (14), (18) and (19) into equations (1) to (9), 
and equating the sums of the coefficients of like powers of r and z to 
zero gives the following system of nonlinear ordinary differential equations: 
du. 


^ = v2-uJ-:^Gpo-rpip-jl'p^-2prr 


dv 


r 

— = - 2u V + ■ v 

dt 15 


^=.w2.±Lcpo-rp,Pi’'p5-2p 


( 20 ) 

( 21 ) 


zz 


(22) 
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dp 


o_ 


dt 

dp 

dp 


= - 2PqUj, - PqW^ 


rr 


— - — “ P(\^r»-rv. “ Pr*r^ 


rr“r '^O’^rrz ^rr’^z 


zz 


dt 

du 


- *2PqUj, 22 " “ ^O^zzz " ^^ZZ^J 


du 


Jrzz ^ 
dt 


= ^v^^'rrr ‘ '*“r"rrr ' ^rir-2)pf-j^pl-^^ 
at 


= 2Vr^rzz * ^r\zz ' ^z'^rzz ' ^^>-’^^iP'iPo'\tPzz 
dv. 


^ = - 4u V „ - 2v u „ 
r rrr r rrr 


dv 


- 2v^u^__ - 2w_v^_^ 
^ r rzz r rzz z rzz 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 


dw 


^ = - 2UfW„z - ^J^TTZ ■ ’'(’'-2)PlP'l’'Po'%rPzz 
dt 


dw 


zzz _ 


dt 


= - *”z^zzz - 3v(y-2)PiPiX' P 


-y«y-3^2 

ZZ 


(31) 


The set of equations (20) to (31) results from retaining terms through 
third order In r and z in equations (11) to (14), (18) and (19). That Is, 
equations (20) to (31) form a third order set of equations. If, on the 
other hand, we retain only terms of first order in r and z, we get 
the following first-order set: 
du 


dt ' 3 ' 


( 32 ) 
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dv 

— ^ = - 2u 
dt 

(33) 

4ff p 

— = - w_ Gp* 

dt * 3 0 

(34) 


(35) 

If we retain terms through second order in eqiiations (11) to (14)^ 

(18), and (16), we get the following second-order set of equations: 

- ,,2 „2 4n __ „-y.y-2. 

-J- - ’'r - "r - 7 - ^Pl^l f>0 »tT 

(36) 

dv 

-ir“-=*Vr 

(37) 

^ . ±L GPq - ypiPl^Pl'^P^ 

(3e 

~ = - 2P(,u, - PgW^ 

(39) 

dP-r 

— = - iP„\ - Prr'^z 

(40) 

dp 

— ^ - - 2p,_u. - 3p__w_ 
*^zz r '^zz z 

(41) 

Several interesting observations can be made by comparison of the first, 
second, and third order sets of equations. Examination of the first order 
set (eqs. (32) to (35)) shows that gravitational and rotational forces (last 
two terms in eq. (32) and last term in (34)) sq>pear when first-order 

terms are retained in the series expansions in r and z. 

The second- 


order set (eqs. (36) to 41)) contains in addition^ pressure gradient terms 
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(last terms in eqs. (36) and (38)). Finally, if third-order terms are 
retained in the series expansions (eqs. (20) to (31)), a tuii)iilent vis- 
cosity term aiqpears (last term in eq, 21)). Another way of saying this 
is that gravitational and rotational terms appear in the equations when 
only first spacial derivatives (u^, v^, and w^) are present at the origin, 
pressure gradient effects require, in addition, second spacial deriva- 
tives (p^^ and and turbulent viscosity effects require the presence 
of the third spacial derivative v^^^ at the origin. 

Before we can solve the set of equations (20) to (31) numerically, 
we must set conditions at an initial time, say at t - 0. We set the initial 
radial and axial flows equal to zero in the vicinity of the origin, so tliat 


0 0 
uj = u: 


= ui„ = = w° . = w“._ = 0 


4. -z 

where the superscripts zero indicate values at t - 0. Also we take the 
angular velocity ^ = as initially uniform near the origin, so that 


= w® - 
rrz ^zzz 


(42) 



(43) 


and 


V? = V? ^ 0 (44) 

rrr rzz 

Finally, we specify the initial density in the vinicity of the origin. 
For doing this, we introduce an initial disturbance of radius r ^ such 
that for z = 0, 

p^ = Ap 



( 45 ) 


and for r = 0, 
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0 A 0 




'i/ j 


+ p. 


(46) 


where Ap^ is the magnitude of the Initial density fluctuation at the origin 
and n is an integer ^2. This initial density distribution is convenient, 
since its shape at the origin can be altered by varying n, and since it 
gives the ambient density p^ at r^. Moreover it gives zero density 
gradient where the density becomes ambient at r^, so that the distribution 
given by equations (45) and (46) joins smoothly with the uniform ambient 
distribution outside of r j . For r = z = 0, 

p° = p° = Ap® + Pj (47) 


and 



in agreement with equation (14). Taking the second derivatives of 
equations (45) and (46) and setting r = z = 0, we get 



for n = 2 


(48) 


and 

pJr ' Piz = 0 *“>■ n > 2 («) 

IV. NUMERICAL RESULTS AND DISCUSSION 
The set of equations (20) to (31) was solved numerically, subject to 
the initial conditions at t = 0 given by equations (42) to (49). Before 
carr]ring out the solution, the equations and initial conditions were con- 
verted to dimensionless form. For example, equation (20) can be 
written as 
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where 



*2 in 


u_- - _ G pq - n>iPo 


♦ *v-2 * 
^ P 


rr 


(50) 


’ u ^ V 

r — , r — 


Uf 


U) 


* 0 * * ^0 
t = W t, Pa = — 


* Gp. 
G : 


Pl = 


2 .2 




Pjrj 


T7P 


l^rr 


rr 


Similar dimensionless equations are obtained from equations (21) to (31) 
and the initial conditions. As before, the superscripts zorc indicate values 
at t s 0, the subscripts zero indicate values at r = z = 0, and the sub- 
scripts one designate constant ambient values. Thus the dimensionless 
dependent variables can be written as functions of a dimensionless time, 
a gravitational parameter, and a pressure parameter. That is, 

u* = f(t*,G*,pJ) (51) 

Similar equations are obtained for the other dependent variables. X>imen- 

sionless quantities other than those in equation (51) can, of course, be 

used in their place, so long as the same total number of variables appears 

* ♦ • 

on the right side. For instance by intermultiplying t , G and p^, we get 
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Gl/2pl/2t = 


CD' 




1 

cn/2 


and 


Pi _Pl 

r«2 2 p* 

Gp^rj G 


$0 that in place of equation (51) we have 




CD 


Pi 


CD 


p 1/2^ 1/2 p^2 2 

G Pi 


(52) 


When written this way, we can think of the dependent variables as fun- 
tions of a dimensionless time, a rotational parameter, and a pressure 
parameter. Note that in equation (52) the rotation is confined to one 
parameter, in contrast to equation (51) where, instead, gravity occurs 
in only one parameter. In both cases the pressure is confined to one 
parameter. 


a) Uniform Initial Density Distribution 
Consider first the case where the initial disturbance at the origin 
Ap^ is zero (see eqs. (45) and (46)). That is, the density is initially 
uniform as in Larson's (1972) rotating case. However, the present case 
differs from that of Larson because he assumed that the radial velocity 
remains zero at an outer boundary. For Ap® = 0, equations (48) and 
(49) indicate p^„ = Pr,„ - 0 initially, regardless of the value of n in 
equations (45) and (46). But if p^^ and p^^, as well as the radial veloc- 
ity are initially zero near the origin (eq. (42)), equations (20) to (31) 
indicate that those quantities will remain zero. Numerical solution of 
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the equations confirmed this; small disturbances (round off errors) did 
not grow with time. Thus, if the initial density is uniform - 0), 
the pressure and turbulent viscosity terms in equations (20) to (31) will 
be zero at all times, and the equations reduce to the first-order set 
given by equations (32) to (35). 

Calculated results for an initially uniform density are plotted in 
figures 2 to 4. A modified Gear method (Windmarsh 1974) was used for 
the numerical computations. In agreement with the preceding discussion 
the results are independent of pressure-gradient and turbulent viscosity 
effects. This is true regardless of the value of the polytropic exponent y. 
The density changes comparatively slowly over a considerable time span 
and then begins to change rapidly. For each value of the rotation param- 
eter there is a particular dimensionless time, designated the collapse 
time, for which the density and other dependent variables at the origin 
increase without limit. That is, the rotating cloud or vortex tends to 
collapse at that time to form a star. The collapse time increases with 
angular velocity w because the centrifugal field produced by the rotation 
tends to prevent collapse. 

Another effect of the rotation is that for large initial angular veloc- 
cities (iPf the density at the center Pq can decrease before it increases, 
as in figure 2. This is again because of the centrifugal field associated 
with the rotation which tends to throw the gas outward, in opposition to 
the gravitational field. Figure 4 shows that for the larger dimensionless 
angular velocity, the radial flow near the center is outward (u^ is posi- 
tive) until shortly before the collapse time. On the other hand the flow 
near the center in the axial direction (z-direction) is always inward 
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(w, 1 a negative) because of gravitational attraction and the absence of 
centrifugal effects in the axial direction. This inward axial flow is in 
fact the main reason that collapse can eventually occur even when the 
radial flow is outward for awhile. When the radial flow is outward, tae 
angular velocity at the origin ca decreases with time as in figure 3. 

This decreases the centrifugal force field so that the radial flow can 
become negative (u^ becomes negative), and collapse eventually occurs. 
But if It were not for the axial inflow during this time, the gravitational 
field would be weakened to such an extent because of the decreased den- 
sity, that collapse would not occur. Although figure 2 indicates that the 
density at the origin Pq can decrease because of the radial outflow, 
that decrease is not nearly as great as is would be if the axial inflow 
were absent. Thus, the collapse process for a rotating cloud can be 
much more complicated that it is for a> = 0. 

In order to investigate the spacial variations of the density and veloc 
tty, and to give somewhat more confidence in the series solution for 
r = z s 0 plotted in figures 2 to 4 an approximate numerical solution of 
the original partial differential equations was obtained. To reduce the 
number of independent spacial variables to one, and thus to hold the re- 
quired computation time within reasonable limits, the solution was ob- 
tained for z = 0. Then w = 0, and we set 


dz Sz r (9u/3r)j.^^.(j 


(53) 


in equation (4), where p had been taken outside of the derivative sign 
because, by symmetry, dp/ dr = 0 at z = 0. Equation (53) is exact if the 
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flow toward or away from the origin is spherically symmetric ((9 w/0z)^_jj_q 
= (9u/3r)j._2_Q). It also gives a consistent result for r = z = 0 for 
spherically nonsymmetric flow, since (equation (11)) u/r approaches 
0u/9r at the origin. Finally it is consistent for no axial flow (0w/0z = 0). 
For (0 w/0z)j._2j_q in equation (53) we use equation (34) since, as men- 
tioned previously, equation (22) reduces to that equation for uniform initial 
density. We take the gravitational potential as spherically symmetric 
in order to give results consistent with those for the series solution. (In 
the series solution <p was taken to be spherically symmetric at the 
origin (eq. (17)). Again we note that Larson’s (1972) work indicates that 
this assumption should give good results. Equation (5) gives, for <p 
spherically symmetric, 




As in the case of the series solution we apply all of the boundary condi- 


tions at the origin. Equations (11) to (14) show that we can use for bound- 

2 2 

ary conditions, u = v = w = 0 v/0r = dp/ 0r = 0, at r = z = 0. 

Equations (1), (2), (4), (6) to (9), (34), (53) and (54) were solved 
numerically by an improved Euler method with increments for r/r j of 
0. 01, except near r = 0 where increments of 0. 005 were used. To give 
results comparable with those for the series solution in figures 2 to 4, 


the density and angular velocity were taken as initially uniform, and u 


as initially zero. (For the series solution those initial conditions were 


applied only in the vicinity of r 
Pl/^Gp^rj) = 0.1, and y = 5/3 


= z = 0). Results for = 1, 

are plotted against r/r^ in figures 5 to 7. 



17 


For comparison with the series solution, densities at r = z = 0 are also 
plotted in figure 2. The agreement is satisfactory. 

Figure 5 shows that the vortex collapses in a nonhomologous (non- 
similar) fashion in agreement with the results of Larson (1972). A 
difference between our results and those of Larson is that the latter showed 
the formation of a ring of mass for rotating flow, rather than a concen- 
tration at the center, (It will be seen later however that our results for 

< 0 are somewhat similar to Larson's results. ) It is not clear 
whether this difference is due to a difference in boundary conditions 
(Larson assumed u = 0 at an outer boundary), to the numerical techniques 
used, or to the inclusion of turbulent viscosity effects (with v = 0 v/9r 
= 0 at the center) in the present calculations. It is of interest that 
Tscharnuter (1975), using conditions similar to those of Larson, but a 
numerical technique differing from both that of Larson and the present 
paper, obtained a concentration of mass at the center. 

The density distributions in figure 5 remain flat near the center of 
the vortex. This flatness is a carry over from the initially uniform 
density distribution. It is in agreement with equations (24), (26), (28) 
and (30), which show that for uniform initial p, u, w, and a>,p will re- 
main uniform near the origin remains zero), as discussed earlier. 
Also figures 6 and 7 show that v and u remain linear near the origin, 
as they should if the first-order set of equations (32) to (35) is to give a 
description of the collapse process. (As discussed previously, the set 
of equations (20) to (31) reduce to the set (32) to (35) in the present case. ) 

Figures 6 and 7 show that the variations of u and v with r are 
nearly linear for a considerable time span (beginning at t = 0). In the 
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case of u this is partly because, for a slowly varying density, the 
gravitational force is nearly linear in r (eq. (54)). The only other term 
in equation (1) which is important at early times is the v /r term, 
which is also nearly linear in r. Thus for early times with p = p^, 
equation (1) integrates to 


u = 




(55) 


and equation (2) gives 




0 

0 ) r 


(56) 


both of which are linear in r. For larger times the variations of u and 

V become highly nonlinear and tend to develop peaks. In the case of u 

this is evidently because the mass becomes concentrated near the origin 

-2 

and the gravitational force is nearly proportional to r outside of the 
mass concentration. Inside of the mass concentration the density is 
uniform so that the gravitational force, and thus u, are still proportional 
to r. The peakedness of the v profiles in the vicinity of the peaks of 
the u curves is due to the fact for large u, v approaches a l/r 
variation (invicid vortex solution). 

To determine the effect of turbulent viscosity on the results, curves 
for that quantity set equal to zero are plotted dashed in figure 6 for com- 
parison with the curves for nonzero turbulent viscosity. The differences 
are slight. In particular the good agreement near the origin indicates 
that the effect of turbulent viscosity on the angular velocity at the origin 
(o) = Vj.) is zero. This is in agreement with the series solution, where 
the turbulent viscosity term drops out for the present case (uniform 
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Initial p). For times larger than those shown, the effect of turbulent 
viscosity in the vicinity of the peak may become greater, but the effect 
on the angular velocity at the origin should still be zero. 

The result that the turbulent viscosity has no effect on the angular 
velocity at the origin may seem to be contrary to experience. For 
instance when the arms of a whirling skater are retracted or extended 
they exert a tangential force on the skater's trunk and thus change the 
angular velocity of the latter. In that case the tangential force is neces* 
sary for changing the angular velocity of the trunk. The difference be- 
tween that case and the rotating cloud (where a tangential force or 
turbulent viscosity is not necessary) appears to be that in the latter the 
radial velocity extends all the way to the center and is zero only at the 
center. In order to check a case which was con^arable to that of the 
skater, u was set equal to zero for r/r^ between 0 and 0.1. It was 
found that the angular velocity at the center changed with time only when 
the turbulent viscosity was nonzero. 

Perhaps the most important effect of turbulent viscosity is that it 
enables the assumed initial solid-body- like rotation to be realized. In the 
absence of a turbulent viscosity the v profile could be arbitrary, and 
there would be no assurance that y and its spacial derivatives are finite 
at the origin, as required for the series solution. The presence of the 
eddy viscosity however provides a tangential stiffness, so that the assumed 
wheel flow can be attained, particularly near the origin. Figure 6 and 
equation (56) show that once an initial wheel flow is established, it can 
remain for a considerable time, even in the presence of small radial 
flows. 
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b) Nonuniform Initial Density Distributions 

Consider first the case where ^ 0 and n > 2 in equations (45) 
and (46). Then equation (49) shows that and p__ are initially zero, 
and equations (20) to (31) show that they will remain zero for the initial 
conditions used herein. Thus, as in the case of Ap^ = 0, the set of third- 
order equations (20) to (31) reduces to the first-order set (32) to (35). 

The results for this case will therefore be similar to those for uniform 
Initial density at the origin, although Pq will be different (for the same 
ambient density p^). The discussion for uniform initial density applies 
to the non-uniform case when n > 2 in equations (45) and (46), at least 
to the present order of approximation. Density evolution curves for 
Ap®/pj = 0. 1 {pq/Pi = 1. ij are plotted in figure 8. As expected, these 
curves, as well as those for the velocity components (not shown), are 
similar to those for Ap^/p^ = 0. They are independent of pressure and 
the polytropic exponent y. 

Dimensionless coll^se times for various values of rotation parameter 
and Ap®/pj are plotted in figure 9. We note that the curves for Ap® * 0 
could be obtained by multiplying the ordinates on the latter by 
^1 + Ap®/pjj'^''^^ and the abscissas by + Ap^/p^ . The collapse 
times are of the same order of magnitude as the free-fall times, in 
agreement with the results of Larson and others. 

To determine the effect of a non-uniform initial density profile on the 
profiles at later times, the governing partial differential equations were 
solved by the method and approximations used for the uniform initial 
density case in figure 5, The results are shown in figure 10, where p/pj 
is plotted against r/r^ for various dimensionless times, and for n = 4 
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and == ±0, 1 in equations (45) and (46). The rest of the parameters 

are the same as Uiose in figure 5. Comparison of figures 5 and 10 indi- 
cates that the initial density profile can have a large effect on the evolution 
of the profiles. In particular the curves for = -0.1 (Fig. 10(b)) 

show the development of a pronounced peak away from r = 0. These 
profiles are somewhat similar to those of Larson (1972) for rotating flow 
and uniform initial density. In the present case the effect does not seem 
to be entirely due to rotation, because when v was set equal to zero, the 
curves, although considerably altered, still showed the development of 
a peak away from r = 0. This concentration of mass in a ring may break 
up into a binary star system (Larson, 1972). For all three cases the values 
of d p/dr at r = 0 remained zero, as they should according to the series 
solution for these cases. 

We still have to consider the case where n = 2, pj > 0, and Ap® * 0 

in equations (45) and (46). This is the only case for which p„^ and p„„, 

rr zz’ 

according to the present third-order solution, are not equal to zero, so 
that we have to consider the full set of equations (20) to (31). The pressure 
and the polytropic exponent y, as well as the turbulent viscosity, may 
have an effect on quantities at r = z = 0 for this case. 

Results for n = 2, ) = 1» ^^/Pj = ±0. 1, and y = 1 and 

5/3 are plotted in figure 11. First, second, and third order approximations 
for the evolution of densities at the origin are shown. (For the previous 
cases where n > 2 or Ap^ = 0, the three approximations were of course 
identical. ) The fact that the second and third order approximations are 
nearly the same, particularly for y = 1 (optically thin case), indicates 
that the third order approximation gives quite accurate results. The 
dashed curves for the first approximation are also the curves for zero 
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pressure (zero temperature for nonzero density) for the three approxi- 
mations. Results for small pressures were essentially the same as 
those for zero pressure, indicating that the results for zero pressure are 
stable. The effect of pressure for positive is to slow up the collapse 
process, since the pressure-gradient force is away from the center. (The 
pressure gradient is in the same direction as the density gradient, accord- 
ing to eq, (6a).) For negative the opposite trend occurs. Increasing 
the value of y or of amplifies these effects. Tne turbulent viscosity 
term in equation (21), a third-order term, was found to be negligibly small 
in most cases. 

The solid curve for y - 5/3 in figure 11(a) indicates that contraction 
stops after a maximum density is reached. The curves for r = 1, as 
well as those for y - 5/3 in figure 11(b), however, show a strong tendency 
for collapse to continue. In the latter cases the gravitational force appar- 
ently increases faster than the resisting (for positive pressure- 
gradient force. These trends, except those for y - 5/3 and a negative 

in figure 11(b), are in agreement with simplified analyses which 
indicate that only for y < 4/3 will collapse occur (e.g. Schatzman 1972). 
For a negative Ap^ the pressure-gradient force can aid, rather than 
hinder the collapse process so that it does not seem surprising that there 
can be a strong tendency for collapse to occur in that case, even for 
y s 5/3 . Of course as discussed earlier, for uniform initial density 
( Ap^ s 0), or for n > 2 in equations (45) and (46), that is, if d^p/dr^ 
at the center is initially zero or negligibly small, our series solution 
indicates that collapse can occur regardless of the values of y and of 
pressure parameter. Thus if the initial disturbance can be approximated 



23 


by equations (45) and (46) the pressure and y, according to the present 
third order analysis, should have an effect on the collapse process at the 
origin only for a rather special case, that is, for n - 2 and ^ * 0. 

I should like to acknowledge the work of Frank Molls on the 
numerical solutions of the differential equations in the present paper. 
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Figura 6. - Davtiopfflaflt ad tmgarrtlal v^ly grofilt In w>rttx. 

9l/IC9frf>-0.l. r* Vi. Density and 
angular velocity m InHiatly unMorm. and radial and ail^ 
velocKlK aro Initially zaro. RtsuHs nav cantor of vortex 
are indapandai^ ot prassure, poiytropic «xponant. ard tur- 
duiwt viscosity aftOcts. 




Figure 7. - Davetepmetd et radial velocity profile in vortex. 
«®/IC‘^^p 1^^IM, P|/ICp{rf).0.1. y-5/3. Density 
and at^lar velocity are initially uniforin, and ridbl 
and aiW vatecHtas are initially nro. Results nMr can- 
tar of «rt«r are indipandint of pressure, potytroiiik m- 
ponant. and tiJiiu lent viscosity affects. 
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(8) • 0.1 in equations (43) and (46). 


Fl^re 10. ' Development of density profile In vortex, n ■ 4 
In equations (45) and (46). u°/(G^'^pp) • 1. PiKGpfr})- 
0.1. y • 5/3i An^lar velocity Is Initially uniform, and 
radial and axial velocities are initially zero. Results near 
center of vortex are independent of pressure, polytropic ex- 
ponent, and turbulent viscosity effects. 



Figure 10. - Concluded. 



01 Aj)®/pi * 0.1 in equations (45) and (46). 

Figure 11. - Time evolution of density at center of vortex for n • 2 in 
equations (45) and (46). Pi/(Gf^rJ) • 1. Angular velocity near cen- 
ter is initially uniform, and radial and axial velocities near center 
are iritially zero. 





